clear D tempROI tempWSS tempti pROI pWSS ti th0 th90 th180 th270

% this program takes the wkspace variables and resamples the ROI image into
% polar space for shear rate/shear stress calculations

mu        = 0.032;  % poise (.032 g/(cm s) for blood, .01 for H20)
method    = SV_parameters.interp;
plotFLAG  = 0;      % show images on the screen
movieFLAG = 0;      % save the images as a movie
phases    = size(VelDironeVelDataMask,3);
V1        = VelDironeVelDataMask(:,:,1);
[M,N]     = size(V1);
limits    = [-N/2 N/2 -M/2 M/2];
pROI      = [0 0];
pWSS      = 0;
ti        = 0;
Th        = (0:360)*pi/180;% this sets the angular limit and resolution of the resample
radialres = 0.5;           % this sets the radial resoltion of the polar resample
if M >= N                  % this sets the radial limit of the polar resample
    radiallimit = .71*M;   
else
    radiallimit = .71*N;
end


if plotFLAG
    if movieFLAG
        hfig = figure('Position',[10 50 400 400],...
            'Color',[1 1 1], 'Name','Cartesian to Polar transformation');
    else
        hfig = figure('Position',[10 50 400 400],...
            'Color',[1 1 1], 'Name','Cartesian to Polar transformation','WindowStyle', 'docked');
    end
    polarimageaxes  = axes('OuterPosition', [0.0121 0.538 0.4318 0.4186]);
    shearstressaxes = axes('OuterPosition', [0.5100 0.5378 0.4318 0.4186]);
    polarstressaxes = axes('OuterPosition', [0.0121 .1 0.4318 0.4186]);
    filteredTzraxes = axes('OuterPosition', [0.5100 .1 0.4318 0.4186]);
end

%% Plot time step 
h = waitbar(0,'Transforming to polar space and calculating WSS');

for i = 1:phases    


    % find non-zero ROI coords:
    x = find(b(:,2,i) > 0);
    y = find(b(:,1,i) > 0);
    roisize = size(x,1);
    centerx=stats(i).Centroid(1);
    centery=stats(i).Centroid(2);
    V1=VelDironeVelDataMask(:,:,i);
    
%     if plotFLAG
%         imageaxes = subplot(2,2,1);
%         imagesc(V1, [-VENC/2 VENC/2]);
%         axis image; hold on
%         plot(centerx, centery,'+')
%         plot(b(y,1,i),b(x,2,i),'w','LineWidth',2);
%         hold off
%     end

%Plot conversion to polar coordinates

    % Calculate polar coordinates of image
    [X Y]=meshgrid((1:N)-centerx,(1:M)-centery);
    Zi(:,:,i) = cart2pol2(X,Y,V1,Th,(0:radialres:radiallimit)',method); % create variable for 38

    % Plot polar re-sample back into cartesian coords
    if plotFLAG
        [Th2,rho2]=meshgrid(Th,(0:radialres:radiallimit)');
        [X,Y]=pol2cart(Th2,rho2);
        axes(polarimageaxes);
        surf(X,Y,Zi(:,:,i),'EdgeColor','none');
        caxis([-VENC/2 VENC/2]);
        view(2); axis square; axis ij;box on; hold on
        plot3(b(y,1,i)-centerx,b(x,2,i)-centery,ones(roisize,1)*VENC,'w','LineWidth',2);
        plot3(0,0,VENC,'+')
        axis(limits)
        colorbar([0.4159 0.5935 0.0192 0.3215],...
            'Box','on',...
            'XLim',[-0.5 1.5],...
            'Location','manual');
        title('Velocity Map, Resampled Polar Coords','FontSize',8)
        xlabel('Pixels','FontSize',8)
        ylabel('Pixels','FontSize',8)
        hold off
    end

% Plot Shear stress in both coordinate systems

    % Calc shear stress in polar space
    [FX FY]    = gradient(Zi(:,:,i),pi/180*structInt.PixelSpacingColumns,radialres*structInt.PixelSpacingRows*.1); % this is incorrect for anglular shear... pixel spacing is in mm need to convert to cm...
    Tzr(:,:,i) = FY.*mu;
    [tempROI tempWSS tempti] = polarROI(Tzr(:,:,i),b(:,:,i),stats(i),radialres);

    % Pad Array if the temp sizes aren't the same as the temporal arrays
    sizetempROI = size(tempROI,1);
    sizetempWSS = size(tempWSS,1);
    sizetempti  = size(tempti,1);
    
    sizediffROI = sizetempROI - size(pROI,1);
    sizediffWSS = sizetempWSS - size(pWSS,1);
    sizediffti  = sizetempti - size(ti,1);
    
    if sizediffROI > 0
        pROI = padarray(pROI,[sizediffROI 0],NaN,'post');
    else
        tempROI = padarray(tempROI,[-sizediffROI 0],NaN,'post');
    end

    if sizediffWSS > 0
        pWSS = padarray(pWSS,[sizediffWSS 0],NaN,'post');
    else
        tempWSS = padarray(tempWSS,[-sizediffWSS 0],NaN,'post');
    end
    
    if sizediffti > 0
        ti = padarray(ti,[sizediffti 0],NaN,'post');
    else
        tempti = padarray(tempti,[-sizediffti 0],NaN,'post');
    end
    
    pROI(:,:,i) = tempROI;
    pWSS(:,i)   = tempWSS;
    ti(:,i)     = tempti;

    % Plot shear stress in polar space   
    if plotFLAG
        axes(polarstressaxes);
        surf((0:360),rho2,Tzr(:,:,i),'EdgeColor','none'); hold on
        caxis([-VENC/10 VENC/10]);
        view(2); axis square; axis xy; box on;
%         axis([0 360 -20 radiallimit])
        axis([0 360 0 radiallimit])
        plot3(pROI(1:sizetempROI,1,i),pROI(1:sizetempROI,2,i)/2+0.5,ones(sizetempROI,1)*VENC/8,'w','Linewidth',2);
        plot3(pROI(1:sizetempROI,1,i),pROI(1:sizetempROI,2,i)/2-0.5,ones(sizetempROI,1)*VENC/8,'w','Linewidth',2);
%         plot3(ti(1:sizetempti,i), pWSS(1:sizetempti,i),ones(sizetempti,1)*VENC/8);
        colorbar1 = colorbar([0.4159 0.152 0.01918 0.3215],...
            'Box','on',...
            'XLim',[-0.5 1.5],...
            'Location','manual');
        set(polarstressaxes,'XTick',[0 90 180 270 360]);
        hold off
        title('WSS Map, Polar Coords','FontSize',8)
        xlabel('Angular Position (\theta)','FontSize',8)
        ylabel('Radial Position (\rho)','FontSize',8)
    end

    % Plot shear stress in cartesian space
    if plotFLAG
        axes(shearstressaxes);
        surf(X,Y,Tzr(:,:,i),'EdgeColor','none')
        caxis([-VENC/10 VENC/10]);
        view(2); axis square; axis ij; box on; hold on
        plot3(0,0,VENC/8,'+')
%         plot3(b(y,1,i)-centerx,b(x,2,i)-centery,ones(roisize,1)*VENC/8,'w','LineWidth',2);
        axis(limits)
        colorbar1 = colorbar([0.915 0.5935 0.01918 0.3215],...
            'Box','on',...
            'XLim',[-0.5 1.5],...
            'Location','manual');
        hold off
        title('WSS Map, Cartesian Coords','FontSize',8)
        xlabel('Pixels','FontSize',8)
        ylabel('Pixels','FontSize',8)
    end
    
    tempWSS = pWSS(~isnan(pWSS(:,i)),i); % remove Nan's
    tempti  = ti(~isnan(pWSS(:,i)),i);
    
    th0   = find(tempti>355 | tempti<5); % find shear @ positions
    th45  = find(tempti>40 & tempti<50);
    th90  = find(tempti>85 & tempti<95);
    th135 = find(tempti>130 & tempti<140);
    th180 = find(tempti>175 & tempti<185);
    th225 = find(tempti>220 & tempti<230);
    th270 = find(tempti>265 & tempti<275);
    th315 = find(tempti>310 & tempti<320);
    
    WSSlocal(1,i) = mean(tempWSS(th0)); 
    WSSlocal(2,i) = mean(tempWSS(th45));
    WSSlocal(3,i) = mean(tempWSS(th90));
    WSSlocal(4,i) = mean(tempWSS(th135));
    WSSlocal(5,i) = mean(tempWSS(th180));
    WSSlocal(6,i) = mean(tempWSS(th225));
    WSSlocal(7,i) = mean(tempWSS(th270));
    WSSlocal(8,i) = mean(tempWSS(th315));
    
    % Plot WSS
    if plotFLAG
        axes(filteredTzraxes);
        plot(ti(1:sizetempti,i), pWSS(1:sizetempti,i),'r');
        hold on;
%         filtered=medfilt1(tempWSS(1:sizetempWSS-1),15);
%         plot(tempti(1:sizetempti-1),filtered,'r')
        axis([0 360 -30 30])
        axis square
        plot(ti([th0;th90;th180;th270],i),pWSS([th0;th90;th180;th270],i),'o');
        set(filteredTzraxes,'XTick',[0 90 180 270 360]);
        grid on;
        hold off;
        title('WSS along ROI line','FontSize',8)
        xlabel('Angular Position (\theta)','FontSize',8)
        ylabel('WSS (dynes/cn^2)','FontSize',8)
    end
    
    % prep for movie
    if movieFLAG && plotFLAG
        tempframe = getframe(hfig);
        D(:,:,:,i) = tempframe.cdata;
    end

    waitbar(i/phases, h);
end
close(h); 
%% Save Movie File

if movieFLAG && plotFLAG
    figure;
    mov = immovie(D);
    movie2avi(mov,'shear.avi','compression','CinePak');
    clear mov D tempframe; close;
end


        

